home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qpsrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.9 KB  |  114 lines
  1. /* integration/qpsrt.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. static inline void 
  21. qpsrt (gsl_integration_workspace * workspace);
  22.  
  23. static inline
  24. void qpsrt (gsl_integration_workspace * workspace)
  25. {
  26.   const size_t last = workspace->size - 1;
  27.   const size_t limit = workspace->limit;
  28.  
  29.   double * elist = workspace->elist;
  30.   size_t * order = workspace->order;
  31.  
  32.   double errmax ;
  33.   double errmin ;
  34.   int i, k, top;
  35.  
  36.   size_t i_nrmax = workspace->nrmax;
  37.   size_t i_maxerr = order[i_nrmax] ;
  38.   
  39.   /* Check whether the list contains more than two error estimates */
  40.  
  41.   if (last < 2) 
  42.     {
  43.       order[0] = 0 ;
  44.       order[1] = 1 ;
  45.       workspace->i = i_maxerr ;
  46.       return ;
  47.     }
  48.  
  49.   errmax = elist[i_maxerr] ;
  50.  
  51.   /* This part of the routine is only executed if, due to a difficult
  52.      integrand, subdivision increased the error estimate. In the normal
  53.      case the insert procedure should start after the nrmax-th largest
  54.      error estimate. */
  55.  
  56.   while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]]) 
  57.     {
  58.       order[i_nrmax] = order[i_nrmax - 1] ;
  59.       i_nrmax-- ;
  60.     } 
  61.  
  62.   /* Compute the number of elements in the list to be maintained in
  63.      descending order. This number depends on the number of
  64.      subdivisions still allowed. */
  65.   
  66.   if(last < (limit/2 + 2)) 
  67.     {
  68.       top = last ;
  69.     }
  70.   else
  71.     {
  72.       top = limit - last + 1;
  73.     }
  74.   
  75.   /* Insert errmax by traversing the list top-down, starting
  76.      comparison from the element elist(order(i_nrmax+1)). */
  77.   
  78.   i = i_nrmax + 1 ;
  79.   
  80.   /* The order of the tests in the following line is important to
  81.      prevent a segmentation fault */
  82.  
  83.   while (i < top && errmax < elist[order[i]])
  84.     {
  85.       order[i-1] = order[i] ;
  86.       i++ ;
  87.     }
  88.   
  89.   order[i-1] = i_maxerr ;
  90.   
  91.   /* Insert errmin by traversing the list bottom-up */
  92.   
  93.   errmin = elist[last] ;
  94.   
  95.   k = top - 1 ;
  96.   
  97.   while (k > i - 2 && errmin >= elist[order[k]])
  98.     {
  99.       order[k+1] = order[k] ;
  100.       k-- ;
  101.     }
  102.   
  103.   order[k+1] = last ;
  104.  
  105.   /* Set i_max and e_max */
  106.  
  107.   i_maxerr = order[i_nrmax] ;
  108.   
  109.   workspace->i = i_maxerr ;
  110.   workspace->nrmax = i_nrmax ;
  111. }
  112.  
  113.  
  114.